// ****************************************************************************
// *     Filename:  DetectorConstruction.cc
// *      Created:  2012-09-14 11:34:20
// *      Version:  1.0
// *     Revision:  0
// *       Author:  Tran Hoai Nam (nam), nam@kuno-g.phys.sci.osaka-u.ac.jp
// * Organization:  Osaka University
// *      License:  Geant4 Software License 
// *                http://geant4.cern.ch/license/LICENSE.html
// *  Description:  
// ****************************************************************************
#include "DetectorConstruction.hh"
#include "MySD.hh"
//Material and geometry
#include "G4Material.hh"
#include "G4NistManager.hh"     
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
#include "G4Sphere.hh"
#include "G4SDManager.hh"
#include "G4Colour.hh"
#include "G4VisAttributes.hh"
#include "G4RotationMatrix.hh"
#include "G4SubtractionSolid.hh"
#include "G4UnionSolid.hh"
#include "G4IntersectionSolid.hh"

#include "G4RunManager.hh"
//Field 
#include "G4FieldManager.hh"
#include "G4TransportationManager.hh"
#include "G4UniformMagField.hh"
#include "G4ChordFinder.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4ClassicalRK4.hh"

#include "globals.hh"

DetectorConstruction::DetectorConstruction() : fpPhysWorld(0)
{
	DefineMaterials();
	SetColors();
	SetDefaultParameters();
	//SetMagField(3*tesla,0,0);
	fpSDMan = G4SDManager::GetSDMpointer();
}

DetectorConstruction::~DetectorConstruction()
{; }

void DetectorConstruction::SetDefaultParameters()
{
	fChamberH = 380*mm;
	fChamberR = 310*mm/2.;
	fChamberT = 5*mm;

	fVPortExL = 50*mm;
	fVPortR = 120/2.*mm; 
	fVPortT = 5*mm;
	fBeamWindowT = 0.1*mm;

	fWorldZ = 3.5*fChamberR;
	fWorldX = 3.5*fChamberR;
	fWorldY = fChamberH + 20*mm;

	fTrigPST = 0.5*mm;
	fTrigPSL = 50*mm;
	fSiL = 50*mm;
	fdESiT = 65*um;
	fESiT = 1500*um;
	fVetoPST = 1*mm;
	fTargetT = 200*um;

	fCollT = 1*mm;

	fGeR = 60*mm/2;
	fGeT = 30*mm;
}

G4VPhysicalVolume* DetectorConstruction::Construct()
{
	// World volume, beam in +z direction
	G4Box* BoxWorld
		= new G4Box("Worldbox",fWorldX/2,fWorldY/2,fWorldZ/2);
	G4LogicalVolume *logWorld;
	logWorld = new G4LogicalVolume(BoxWorld, fpVacuum,"Worldlog",0,0,0); 
	fpPhysWorld = new G4PVPlacement(0,G4ThreeVector(),
			logWorld,"expHall",0,false,0);
	logWorld->SetVisAttributes(G4VisAttributes::Invisible); //Set invisible
	
	// Vacuum chamber
	BuildChamber();
	G4RotationMatrix* rm2 = new G4RotationMatrix();
  rm2->rotateX(90*deg);
	new G4PVPlacement(rm2,G4ThreeVector(0.,0.,0.),fpLogChamber,
			"physChamber",logWorld,false,0);
	//fpLogChamber->SetVisAttributes(G4VisAttributes::Invisible);

	G4double posBWindow = -fChamberR - fVPortExL - fBeamWindowT/2;
	new G4PVPlacement (0, G4ThreeVector(0,0,posBWindow),
			fpLogBeamWindow, "physBWindow",logWorld,false,0);

	// Detectors and target
  BuildDetectors();

	G4double posTrig = posBWindow - fBeamWindowT/2 - fTrigPST/2 - 5*mm;
	new G4PVPlacement(0, G4ThreeVector(0,0,posTrig), 
			fpLogTrigPS, "physTrigPS",logWorld,false,0);

	G4double posTrig2 = fChamberR - fChamberT - 10*mm;
	new G4PVPlacement(0, G4ThreeVector(0,0,posTrig2), 
			fpLogTrigPS, "physTrigPS",logWorld,false,1);

	G4RotationMatrix *rm3 = new G4RotationMatrix(); 
	rm3->rotateY(-45*deg);
	new G4PVPlacement(rm3, G4ThreeVector(0,0,0),
			fpLogTarget, "physTarget",logWorld,false,0);

	// Det block
	G4double detBlockL = fSiL + 1*mm;
	G4double gap = 1*mm;
	G4double detBlockT = fVetoPST + fdESiT + fESiT + 3*gap;
	//G4double posDetBlock = (fChamberR - 20*mm)*1.4142/2;
	G4double posDetBlock = 120*mm*1.4142/2;

	G4Box *boxDetBlock = 
		new G4Box("boxDetBlock",detBlockL/2,detBlockL/2, detBlockT/2);
	G4LogicalVolume *logDetBlock0 = 
		new G4LogicalVolume(boxDetBlock,fpVacuum,"logDetBlock0");
	logDetBlock0->SetVisAttributes(G4VisAttributes::Invisible); //Set invisible
	new G4PVPlacement(rm3,G4ThreeVector(-posDetBlock,0,-posDetBlock),
			logDetBlock0,"physDetBlock0",logWorld,false,0);

	G4double posVetoPS1 = -detBlockT/2 + fVetoPST/2 +gap/2;
	new G4PVPlacement(0,G4ThreeVector(0,0,posVetoPS1),
			fpLogVetoPS,"physVetoPS1",logDetBlock0,false,0);

	G4double posESi1 = posVetoPS1 + fVetoPST/2 + gap + fESiT/2;
	new G4PVPlacement(0,G4ThreeVector(0,0,posESi1),
			fpLogESi,"physESi1",logDetBlock0,false,0);

	G4double posdESi1 = posESi1 + fESiT/2 + gap + fdESiT/2;
	new G4PVPlacement(0,G4ThreeVector(0,0,posdESi1),
			fpLogdESi,"physdESi1",logDetBlock0,false,0);

	// Second det block
	G4LogicalVolume *logDetBlock1 = 
		new G4LogicalVolume(boxDetBlock,fpVacuum,"logDetBlock1");
	logDetBlock1->SetVisAttributes(G4VisAttributes::Invisible); //Set invisible
	new G4PVPlacement(rm3,G4ThreeVector(posDetBlock,0,posDetBlock),
			logDetBlock1,"physDetBlock1",logWorld,false,0);

	posVetoPS1 = detBlockT/2 - fVetoPST/2 - gap/2;
	new G4PVPlacement(0,G4ThreeVector(0,0,posVetoPS1),
			fpLogVetoPS,"physVetoPS1",logDetBlock1,false,1);

	posESi1 = posVetoPS1 - fVetoPST/2 - gap - fESiT/2;
	new G4PVPlacement(0,G4ThreeVector(0,0,posESi1),
			fpLogESi,"physESi1",logDetBlock1,false,1);

	posdESi1 = posESi1 - fESiT/2 - gap - fdESiT/2;
	new G4PVPlacement(0,G4ThreeVector(0,0,posdESi1),
			fpLogdESi,"physdESi1",logDetBlock1,false,1);

	// Ge 
	G4RotationMatrix *rm4 = new G4RotationMatrix();
	rm4->rotateY(90*deg);
	new G4PVPlacement(rm4,G4ThreeVector(-fChamberR - fVPortExL - 20*mm,0,0),
			fpLogGe,"physGe",logWorld,false,0);

	// Collimator
	//new G4PVPlacement(0,G4ThreeVector(-70*mm,0,-120*mm),
			//fpLogColl,"physColl",logWorld,false,0);

	//Sensitive detector definition
	AddSensitiveDetectors();

	// Dummy box, for magnetic field setting
	//G4Box* BoxDummy = new G4Box("DummyBox",fDet1X/2, fDet1Y/2, fDet1Z*5);
	//G4LogicalVolume* logDummy = new G4LogicalVolume(BoxDummy,fpVacuum,"logDummy");
	//logDummy->SetFieldManager(fpFieldMan,false); 
	//true: force field to all daughters
	//new G4PVPlacement(0,G4ThreeVector(0.,0.,0.),
			//logDummy,"physDummy",logWorld,false,0);
  //logDummy->SetVisAttributes(fpYellow);

	// Finish
	return fpPhysWorld;
}

void DetectorConstruction::AddSensitiveDetectors()
{
	G4VSensitiveDetector* trigPSSD = new MySD("trigPSSD"); 
	fpSDMan->AddNewDetector(trigPSSD);
	fpLogTrigPS->SetSensitiveDetector(trigPSSD);

	G4VSensitiveDetector* targetSD = new MySD("targetSD"); 
	fpSDMan->AddNewDetector(targetSD);
	fpLogTarget->SetSensitiveDetector(targetSD);

	G4VSensitiveDetector* vetoPSSD = new MySD("vetoPSSD"); 
	fpSDMan->AddNewDetector(vetoPSSD);
	fpLogVetoPS->SetSensitiveDetector(vetoPSSD);

	G4VSensitiveDetector* chamberSD = new MySD("chamberSD"); 
	fpSDMan->AddNewDetector(chamberSD);
	fpLogChamber->SetSensitiveDetector(chamberSD);

	G4VSensitiveDetector* eSiSD = new MySD("eSiSD"); 
	fpSDMan->AddNewDetector(eSiSD);
	fpLogESi->SetSensitiveDetector(eSiSD);

	G4VSensitiveDetector* dESiSD = new MySD("dESiSD"); 
	fpSDMan->AddNewDetector(dESiSD);
	fpLogdESi->SetSensitiveDetector(dESiSD);
}

void DetectorConstruction::DefineMaterials()
{
	G4NistManager* man = G4NistManager::Instance();
	fpVacuum = man->FindOrBuildMaterial("G4_Galactic");
	fpAir = man->FindOrBuildMaterial("G4_AIR");
	fpCopper = man->FindOrBuildMaterial("G4_Cu");
	fpAluminum = man->FindOrBuildMaterial("G4_Al");
	fpNaI = man->FindOrBuildMaterial("G4_SODIUM_IODIDE");
	fpSilicon = man->FindOrBuildMaterial("G4_Si");
	fpLead = man->FindOrBuildMaterial("G4_Pb");
	fpPScinti = man->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
	fpMylar = man->FindOrBuildMaterial("G4_MYLAR");
	fpGermanium = man->FindOrBuildMaterial("G4_Ge");

	G4double density;
	G4int ncomponents;
	G4double fractionmass;
	fpSSteel = new G4Material("StainlessSteel", 
			density= 8.06*g/cm3, ncomponents=6);
	G4Element *C = man->FindOrBuildElement("C");
	G4Element *Si = man->FindOrBuildElement("Si");
	G4Element *Cr = man->FindOrBuildElement("Cr");
	G4Element *Mn = man->FindOrBuildElement("Mn");
	G4Element *Fe = man->FindOrBuildElement("Fe");
	G4Element *Ni = man->FindOrBuildElement("Ni");
	fpSSteel->AddElement(C, fractionmass=0.001);
	fpSSteel->AddElement(Si, fractionmass=0.007);
	fpSSteel->AddElement(Cr, fractionmass=0.18);
	fpSSteel->AddElement(Mn, fractionmass=0.01);
	fpSSteel->AddElement(Fe, fractionmass=0.712);
	fpSSteel->AddElement(Ni, fractionmass=0.09);
}

void DetectorConstruction::BuildChamber() 
{
	// Vacuum chamber
	G4Tubs *tubBPort = new G4Tubs("tubBPort",
			0,fVPortR+fVPortT,fVPortExL,0*deg,360*deg);
	G4Tubs* tubChamber = new G4Tubs("tubChamber",
			0, fChamberR+fChamberT, fChamberH/2, 0*deg, 360*deg);
	G4Tubs *tubCut1 = new G4Tubs("tubCut1",
			0,fChamberR,fChamberH/2 - fChamberT,0*deg,360*deg);
	G4Tubs *tubCut2 = new G4Tubs("tubCut2",
			0,fVPortR,fVPortExL + 1*cm,0*deg,360*deg);

  G4RotationMatrix* rm = new G4RotationMatrix();
  rm->rotateX(-90*deg);

	G4VSolid *boolSolid1 = new G4UnionSolid("boolSolid1",tubChamber,tubBPort,
			rm,G4ThreeVector(0,fChamberR,0));
	G4VSolid *boolSolid2 = new G4SubtractionSolid("boolSolid2",boolSolid1,tubCut2,
			rm,G4ThreeVector(0,fChamberR,0));

  //rm->rotateY(45*deg);
	//G4VSolid *boolSolid3 = new G4UnionSolid("boolSolid3",boolSolid2,tubBPort,
			//rm,G4ThreeVector(-fChamberR*1.4142/2,-fChamberR*1.4142/2,0));
	//G4VSolid *boolSolid4 = new G4SubtractionSolid("boolSolid4",boolSolid3,tubCut2,
			//rm,G4ThreeVector(-fChamberR*1.4142/2,-fChamberR*1.4142/2,0));

	rm->rotateY(90*deg);
	G4VSolid *boolSolid3 = new G4UnionSolid("boolSolid3",boolSolid2,tubBPort,
			rm,G4ThreeVector(-fChamberR,0,0));
	G4VSolid *boolSolid4 = new G4SubtractionSolid("boolSolid4",boolSolid3,tubCut2,
			rm,G4ThreeVector(-fChamberR,0,0));
	G4VSolid *chamberSolid = new G4SubtractionSolid("chamber",boolSolid4,tubCut1,
			0,G4ThreeVector(0,0,0));

	fpLogChamber = new G4LogicalVolume(chamberSolid, fpSSteel,"logChamber");
	fpLogChamber->SetVisAttributes(fpGrey);

	// Beam window
	G4Tubs *tubBWindow = new G4Tubs("tubBWindow",
			0,fVPortR+fVPortT,fBeamWindowT/2,0*deg,360*deg);
	fpLogBeamWindow = new G4LogicalVolume(tubBWindow,fpMylar,"logBWindow");
	fpLogBeamWindow->SetVisAttributes(fpMagenta);
}

void DetectorConstruction::BuildDetectors() 
{
	G4Box *boxTrigPS = new G4Box("boxTrigPS",fTrigPSL/2,fTrigPSL/2,fTrigPST/2);
	fpLogTrigPS = new G4LogicalVolume(boxTrigPS,fpPScinti,"logTrigPS");
	fpLogTrigPS->SetVisAttributes(fpGreen);

	G4Box *boxVetoPS = new G4Box("boxVetoPS",fSiL/2,fSiL/2,fVetoPST/2);
	fpLogVetoPS = new G4LogicalVolume(boxVetoPS,fpPScinti,"logVetoPS");
	fpLogVetoPS->SetVisAttributes(fpCyan);

	G4Box *boxESi = new G4Box("boxESi",fSiL/2,fSiL/2,fESiT/2);
	fpLogESi = new G4LogicalVolume(boxESi,fpSilicon,"logESi");
	fpLogESi->SetVisAttributes(fpCyan);

	G4Box *boxdESi = new G4Box("boxdESi",fSiL/2,fSiL/2,fdESiT/2);
	fpLogdESi = new G4LogicalVolume(boxdESi,fpSilicon,"logdESi");
	fpLogdESi->SetVisAttributes(fpRed);

	G4Box *boxTarget = new G4Box("boxTarget",fSiL/2,fSiL/2,fTargetT/2);
	fpLogTarget = new G4LogicalVolume(boxTarget,fpSilicon,"logTarget");
	fpLogTarget->SetVisAttributes(fpRed);

	G4Box *boxColl = new G4Box("boxColl",fSiL/2,fSiL/2,fCollT/2);
	fpLogColl = new G4LogicalVolume(boxColl,fpLead,"logColl");
	fpLogColl->SetVisAttributes(fpMagenta);

	G4Tubs *tubGe = new G4Tubs("tubGe",0,fGeR,fGeT/2,0*deg,360*deg);
	fpLogGe = new G4LogicalVolume(tubGe,fpGermanium,"logGe");
	fpLogGe->SetVisAttributes(fpMagenta);
}

void DetectorConstruction::SetMagField(G4double Bx, G4double By, G4double Bz)
{
	G4UniformMagField* myField =
		new G4UniformMagField(G4ThreeVector(Bx,By,Bz));//Bx
	G4Mag_UsualEqRhs* myEquation = new G4Mag_UsualEqRhs(myField);
	G4MagIntegratorStepper* myStepper = new G4ClassicalRK4(myEquation);
	G4ChordFinder* myChordFinder = new G4ChordFinder(myField,0.01*mm,myStepper);
	fpFieldMan = new G4FieldManager(myField,myChordFinder);
}

void DetectorConstruction::SetColors() 
{
	fpYellow = new G4VisAttributes(G4Colour::Yellow());
	fpGreen = new G4VisAttributes(G4Colour::Green());
	fpRed = new G4VisAttributes(G4Colour::Red());
	fpCyan = new G4VisAttributes(G4Colour::Cyan());
	fpMagenta = new G4VisAttributes(G4Colour::Magenta());
	fpGrey = new G4VisAttributes(G4Colour::Grey());
}

